<!DOCTYPE html>
<html>
<head><meta charset="utf-8" />
<title>StokesisoP2P0femrate</title><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>

<style type="text/css">
/* Overrides of notebook CSS for static HTML export */
body {
  overflow: visible;
  padding: 8px;
}
div#notebook {
  overflow: visible;
  border-top: none;
}@media print {
  div.cell {
    display: block;
    page-break-inside: avoid;
  } 
  div.output_wrapper { 
    display: block;
    page-break-inside: avoid; 
  }
  div.output { 
    display: block;
    page-break-inside: avoid; 
  }
}
</style>

<!-- Custom stylesheet, it must be in the parent directory as the html file -->
<link rel="stylesheet" type="text/css" media="all" href="../doc.css" />
<link rel="stylesheet" type="text/css" media="all" href="doc.css" />

<!-- Loading mathjax macro -->
<!-- Load mathjax -->
    <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS_HTML"></script>
    <!-- MathJax configuration -->
    <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
        tex2jax: {
            inlineMath: [ ['$','$'], ["\\(","\\)"] ],
            displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
            processEscapes: true,
            processEnvironments: true
        },
        // Center justify equations in code and markdown cells. Elsewhere
        // we use CSS to left justify single line equations in code cells.
        displayAlign: 'center',
        "HTML-CSS": {
            styles: {'.MathJax_Display': {"margin": 0}},
            linebreaks: { automatic: true }
        }
    });
    </script>
    <!-- End of mathjax configuration --></head>
<body>
  <div tabindex="-1" id="notebook" class="border-box-sizing">
    <div class="container" id="notebook-container">

<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="isoP2P0-Finite-Elements-for-Stokes-Equations">isoP2P0 Finite Elements for Stokes Equations<a class="anchor-link" href="#isoP2P0-Finite-Elements-for-Stokes-Equations">&#182;</a></h1>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>This example is to show the convergence of isoP2-P0 finite elements for the Stokes equation on the unit square:</p>
$$- \Delta u + {\rm grad}\, p  = f \quad {\rm div}\, u    = 0  \quad  \text{ in } \quad \Omega,$$<p>with the pure Dirichlet boundary condition. The solver is based on a DGS type smoother.</p>
<p><strong>References</strong>:</p>
<ul>
<li><a href="http://www.math.uci.edu/~chenlong/226/FEMStokes.pdf">Finite Element Methods For Stokes Equations</a></li>
<li><a href="../project/projectFEM.html">Project: Finite Element Methods for Stokes Equations</a></li>
</ul>
<p><strong>Subroutines</strong>:</p>

<pre><code>- StokesisoP2P0
- squareStokes
- femStokes
- Stokesfemrate

</code></pre>
<p>The method is implemented in <code>StokesisoP2P0</code> subroutine and can be tested in <code>squareStokes</code>. Together with other elements (P2P0, P2P1, isoP2P0, isoP2P1, P1bP1), <code>femStokes</code> provides a concise interface to solve Stokes equation. The isoP2-P0 element is tested in <code>Stokesfemrate</code>.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="isoP2-P0-element">isoP2-P0 element<a class="anchor-link" href="#isoP2-P0-element">&#182;</a></h2><p>The velocity is P1 Lagrange element but on a uniform refined mesh and the pressure is P0 piecewise constant element. Therefore this element is very easy to implement.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[3]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">imatlab_export_fig</span><span class="p">(</span><span class="s">&#39;print-png&#39;</span><span class="p">)</span>  <span class="c">% Static png figures.</span>
<span class="n">node</span> <span class="p">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mf">0.5</span><span class="p">,</span> <span class="nb">sqrt</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">/</span><span class="mi">2</span><span class="p">];</span>
<span class="n">elem</span> <span class="p">=</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">];</span>
<span class="n">set</span><span class="p">(</span><span class="n">gcf</span><span class="p">,</span><span class="s">&#39;Units&#39;</span><span class="p">,</span><span class="s">&#39;normal&#39;</span><span class="p">);</span> 
<span class="n">set</span><span class="p">(</span><span class="n">gcf</span><span class="p">,</span><span class="s">&#39;Position&#39;</span><span class="p">,[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mf">0.3</span><span class="p">,</span><span class="mf">0.2</span><span class="p">]);</span>
<span class="n">subplot</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">);</span> <span class="n">showmesh</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span> <span class="n">findelem</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">uniformrefine</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">subplot</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">);</span> <span class="n">showmesh</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span> <span class="n">findnode</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="s">&#39;all&#39;</span><span class="p">,</span><span class="s">&#39;noindex&#39;</span><span class="p">);</span>
</pre></div>

    </div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Dirichlet-boundary-condition">Dirichlet boundary condition<a class="anchor-link" href="#Dirichlet-boundary-condition">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[4]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">%% Setting</span>
<span class="c">% mesh</span>
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">squaremesh</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span><span class="mf">0.25</span><span class="p">);</span>
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">uniformrefine</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">bdFlag</span> <span class="p">=</span> <span class="n">setboundary</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;Dirichlet&#39;</span><span class="p">);</span>
<span class="n">mesh</span> <span class="p">=</span> <span class="n">struct</span><span class="p">(</span><span class="s">&#39;node&#39;</span><span class="p">,</span><span class="n">node</span><span class="p">,</span><span class="s">&#39;elem&#39;</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;bdFlag&#39;</span><span class="p">,</span><span class="n">bdFlag</span><span class="p">);</span>
<span class="c">% pde</span>
<span class="n">pde</span> <span class="p">=</span> <span class="n">Stokesdata1</span><span class="p">;</span> 
<span class="c">% options</span>
<span class="n">option</span><span class="p">.</span><span class="n">L0</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">maxIt</span> <span class="p">=</span> <span class="mi">4</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">printlevel</span> <span class="p">=</span> <span class="mi">1</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">solver</span> <span class="p">=</span> <span class="s">&#39;mg&#39;</span><span class="p">;</span>
</pre></div>

    </div>
</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[5]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">option</span><span class="p">.</span><span class="n">elemType</span> <span class="p">=</span> <span class="s">&#39;isoP2P0&#39;</span><span class="p">;</span>
<span class="n">femStokes</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">option</span><span class="p">);</span>
</pre></div>

    </div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

    <div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>#dof:   2434,  #nnz:  16930, level:  3  MG WCYCLE iter: 13,  err = 4.5294e-09,  time = 0.13 s
#dof:   9986,  #nnz:  70690, level:  4  MG WCYCLE iter: 13,  err = 8.1896e-09,  time = 0.22 s
#dof:  40450,  #nnz: 288802, level:  5  MG WCYCLE iter: 14,  err = 3.0308e-09,  time =  1.1 s
#dof: 162818,  #nnz: 1167394, level:  6  MG WCYCLE iter: 14,  err = 3.4657e-09,  time =  4.9 s
Table: Error
 #Dof        h      |u_I-u_h|_1    ||u-u_h||   ||u_I-u_h||_{max}

  2690   6.25e-02   6.25359e-01   1.17716e-02   3.30491e-02
 10498   3.12e-02   3.19111e-01   3.05876e-03   8.93540e-03
 41474   1.56e-02   1.61122e-01   7.78951e-04   2.34901e-03
164866   7.81e-03   8.09430e-02   1.96476e-04   6.08646e-04

 #Dof        h      ||p_I-p_h||    ||p-p_h||   

  2690   6.25e-02   2.95825e-01   7.56862e-01
 10498   3.12e-02   1.43194e-01   3.76809e-01
 41474   1.56e-02   7.05586e-02   1.88036e-01
164866   7.81e-03   3.50719e-02   9.39433e-02

Table: CPU time
 #Dof    Assemble     Solve      Error      Mesh    

  2690   7.00e-02   1.30e-01   4.00e-02   1.00e-02
 10498   3.00e-02   2.20e-01   3.00e-02   0.00e+00
 41474   9.00e-02   1.08e+00   1.00e-01   0.00e+00
164866   3.50e-01   4.87e+00   1.70e-01   2.00e-02

</pre>
</div>
</div>

<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Conclusion">Conclusion<a class="anchor-link" href="#Conclusion">&#182;</a></h2><p>Optimal order convergence of velocity and pressure is observed. First order for velocity in H1 norm and for pressure in L2 norm. Second order for velocity in L2 and maximum norm.</p>
<p>Multigrid solver based on DGS smoother converges uniformly.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>To-Do: test other boundary conditions.</p>

</div>
</div>
</div>
    </div>
  </div>
</body>

 


</html>
